function [FEVD, FEVD_opt] = VARfevd(VAR,nsteps,ident)
% =======================================================================
% Compute FEVDs for a VAR model estimated with VARmodel function.
% =======================================================================
% [FEVD, FEVD_opt] = VARfevd(VAR,nsteps,ident)
% -----------------------------------------------------------------------
% INPUT
%   VAR     : structure, result of VARmodel function
%   nsteps  : number of nsteps to compute the IRFs
%
% OPTIONAL INPUT
%   ident   : type of identification: 'oir' (default) or 'bq'
%
% OUTPUT
%   FEVD(t,j,k) : matrix with 't' steps, the FEVD due to 'j' shock for 'k' variable
%   FEVD_opt
%     - FEVD_opt.nsteps : number of steps
%     - FEVD_opt.ident  : identification chosen
%     - FEVD_opt.SE(t,j): matrix with 't' steps, standard error of the 'j' variable
%     - FEVD_opt.invA   : structural matrix such that invA*epsilon = u
% =======================================================================
% Ambrogio Cesa Bianchi, July 2012
% ambrogio.cesabianchi@gmail.com

% I thank Dora Xia for pointing out a typo in the above description.

%% Check inputs
%===============================================
if ~exist('ident','var')
    ident = 'oir';
end

%% Retrieve parameters and preallocate variables
%===============================================
c_case   = VAR.c_case;
nvars    = VAR.neqs;
nlags    = VAR.nlag;
beta     = VAR.beta;
sigma    = VAR.sigma;

MSE   = zeros(nvars,nvars,nsteps);
MSE_j = zeros(nvars,nvars,nsteps);
PSI   = zeros(nvars,nvars,nsteps);
FEVD  = zeros(nsteps,nvars,nvars);
SE    = zeros(nsteps,nvars);

%% Compute the multipliers
%=========================

VARjunk = VAR;
VARjunk.sigma = eye(nvars);

IRFjunk = VARir(VARjunk,nsteps);

% this loop is to save the multipliers for each period
for mm = 1:nvars
    PSI(:,mm,:) = reshape(IRFjunk(:,:,mm)',1,nvars,nsteps);
end

%% Compute the companion matrix (just for BQ identification)
%===========================================================
switch c_case
    case 0
        beta_t = beta';
        F = [beta_t(:,1:nvars*nlags)
            eye(nvars*(nlags-1)) zeros(nvars*(nlags-1),nvars)];
    case 1
        beta_t = beta';
        F = [beta_t(:,2:nvars*nlags+1)
            eye(nvars*(nlags-1)) zeros(nvars*(nlags-1),nvars)];
     case 2
        beta_t = beta';
        F = [beta_t(:,3:nvars*nlags+2)
            eye(nvars*(nlags-1)) zeros(nvars*(nlags-1),nvars)];
end

% Compute matrices for BQ identification
Finf_big = inv(eye(length(F))-F); % from the companion
Finf     = Finf_big(1:nvars,1:nvars);
D        = chol(Finf*VAR.sigma*Finf')'; 


%% Calculate the contribution to the MSE for each shock (i.e, FEVD)
%==================================================================

for mm = 1:nvars % loop for the shocks
    
    % Calculate Total Mean Squared Error
    MSE(:,:,1) = sigma;

    for kk = 2:nsteps;
       MSE(:,:,kk) = MSE(:,:,kk-1) + PSI(:,:,kk)*sigma*PSI(:,:,kk)';
    end;
    
    % Get the matrix inv(A) containing the structural impulses depending on
    % the identification scheme chosen
    if strcmp(ident,'oir')
        invA = chol(sigma)';
    elseif strcmp(ident,'bq')
        invA = Finf\D;
    end
    
    % Get the column of invA corresponding to the mm_th shock
    column = invA(:,mm);
    
    % Compute the mean square error
    MSE_j(:,:,1) = column*column';
    for kk = 2:nsteps
        MSE_j(:,:,kk) = MSE_j(:,:,kk-1) + PSI(:,:,kk)*(column*column')*PSI(:,:,kk)';   
    end

    % Compute the Forecast Error Covariance Decomposition
    FECD = MSE_j./MSE;
    
    % Select only the variance terms
    for nn = 1:nsteps
        for ii = 1:nvars
            FEVD(nn,mm,ii) = FECD(ii,ii,nn);
            SE(nn,:) = sqrt( diag(MSE(:,:,nn))' );
        end
    end
end
 
FEVD_opt.nsteps = nsteps;
FEVD_opt.ident  = ident;
FEVD_opt.SE     = SE;
FEVD_opt.invA   = invA;